Resumen

En este trabajo, se describe una aplicación práctica de las Máquinas de Soporte Vectorial, como soporte para acompañar el desarrollo teórico realizado en el primer video de la serie. Se exploran alternativas de funciones kernel, impacto de los parámetros y rendimiento de clasificadores en el caso de estudio: identificación prematura de pacientes con enfermedad PDAC. Se comparan los algoritmos con un modelo de regresión logstica.

Introducción

El cáncer pancreático es una neoplasia extremadamente mortal, con una incidencia del 30% de los tumores diagnosticados a nivel global durante el año 2020. Una vez diagnosticado, la tasa de supervivencia en los siguientes 5 años es de apenas el 10%. No obstante, si el mismo es detectado en etapas tempranas, la probabilidad de supervivencia aumenta. Desafortunadamente, muchos casos de cáncer pancreático no muestran síntomas detectables sino hasta que el mismo ha generado metástasis. En este sentido, un diagnóstico temprano de los pacientes con cáncer pancreático representa una necesidad clínica insatisfecha y una ventana de oportunidad para el desarrollo de nuevas herramientas que asistan a la identificación temprana de esta afección.

Por tal motivo, se eligió aplicar una clasificación supervisada con Máquinas de Soporte vectorial sobre un conjunto de datos publicado recientemente (Debernardi y col, 2020), con el objeto de ilustrar la fortaleza de dicho algoritmo para clasificar pacientes en función de diversos marcadores urinarios, previo al diagnóstico de PDAC.

Base de datos

La base de datos de Debernardi y colaboradores (doi: 10.1371/journal.pmed.1003489.) cuenta con la edad (años), sexo, valores de los biomarcadores urinarios creatinina [mg/ml], LYVE1 [ng/ml], REG1B [ng/ml], y TFF1 [ng/ml] para pacientes en tres categorías de la variable diagnóstico: sanos, con alguna condición benigna, y con PDAC (adenocarcinoma pancreático ductal, el tipo de cáncer pancreático más recurrente). Los datos fueron generados en 4 centros de investigación. En este trabajo se decidió utilizar los datos de pacientes cuyo diagnóstico era “normal” (sanos) o “maligno” (PDAC), conformando una base de datos de 381 registros, y se utilizó la totalidad de los datos sin tener en cuenta su proveniencia.

Fuente original de los datos: https://doi.org/10.1371/journal.pmed.1003489.s009

Procesamiento de base de datos

Carga de librerías

knitr::opts_chunk$set(echo = TRUE)

library(mlr)
## Loading required package: ParamHelpers
## Warning message: 'mlr' is in 'maintenance-only' mode since July 2019.
## Future development will only happen in 'mlr3'
## (<https://mlr3.mlr-org.com>). Due to the focus on 'mlr3' there might be
## uncaught bugs meanwhile in {mlr} - please consider switching.
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
#remotes::install_github("vqv/ggbiplot")
library(ggbiplot)
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## Loading required package: scales
## Loading required package: grid
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(ggpubr)
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
library(knitr)
library(googlesheets)
library(tidyr)
library(plotly)
## 
## Attaching package: 'plotly'
## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Creación de tema general para los gráficos

theme <- theme(text = element_text(size=10),plot.title = element_text(size=12, face="bold.italic",
               hjust = 0.5), axis.title.x = element_text(size=10, face="bold", colour='black'),
               axis.title.y = element_text(size=10, face="bold"),panel.border = element_blank(),
               panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.title = element_text(face="bold"))

Carga de datos

id <- '11wXwKESJrMw6_hqljbOzpn7ncXPjNWeV' # nombre del archivo en google drive
data <- read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download", id))
#data <- read.csv('pancreas.csv')

Preprocesamiento

# Cambio el nombre de las variables a español
colnames(data) <- c('id','cohorte','origen','edad','sexo','diagnosis','estadio','diagnosis_benigno','plasma_CA19_9','creatinina','LYVE1','REG1B','TFF1','REG1A')
# Elijo normal y maligno, y renombro 1=normal, 3=maligno
data <- data %>% filter ( diagnosis=='1' | diagnosis=='3')
data <- data %>% mutate(diagnosis = sub("1", "normal", diagnosis))
data <- data %>% mutate(diagnosis = sub("3", "maligno", diagnosis))
# Transformo como factor
data$sexo <- as.factor(data$sexo)
data$diagnosis <- factor(data$diagnosis, levels= c('normal','maligno'))
# Me quedo sólo con columnas: edad,sexo,diagnosis,creatinina,LYVE1,REG1B,TFF1 y las reordeno
data <- data %>% dplyr::select(6,5,4,10:13) 
# Normalizo variables 
#data <- data %>% mutate_at(c('edad', 'creatinina','LYVE1','REG1B','TFF1'), ~(scale(.) %>% as.vector))
# Imprimo tabla
kable(data.frame(variable = names(data),
           class = sapply(data, class),
           primeros_valores = sapply(data, function(x) paste0(head(x),  collapse = ", ")),
           row.names = NULL))
variable class primeros_valores
diagnosis factor normal, normal, normal, normal, normal, normal
sexo factor F, F, M, M, M, M
edad integer 33, 81, 51, 61, 62, 53
creatinina numeric 1.83222, 0.97266, 0.78039, 0.70122, 0.21489, 0.84825
LYVE1 numeric 0.8932192, 2.037585, 0.1455889, 0.00280488, 0.00085956, 0.003393
REG1B numeric 52.94884, 94.46703, 102.366, 60.579, 65.54, 62.126
TFF1 numeric 654.282174, 209.48825, 461.141, 142.95, 41.088, 59.793

Análisis exploratorio de las variables

Análisis exploratorio de las variables y correlaciones de acuerdo al valor de la variable diagnóstico (gráficos)

# Gráfico ggpairs
data%>% ggpairs(.,mapping=ggplot2::aes(color = diagnosis,alpha = 0.1),
        upper = list(continuous = wrap("cor", size = 2.5),discrete = "blank", combo="blank"),
        lower = list(combo = "box"),progress = F)+
        theme+
        labs(title= 'Descripción de variables en la base de datos', x='Variable', y='Variable')+
        scale_fill_manual(values=c('royalblue2','red'))+
        scale_color_manual(values=c('royalblue2','#ff7474ff'))

Resultados

Separación de conjunto de entrenamiento (train) y prueba (test) [70:30]

set.seed(1409) # para asegurar reproducibilidad
dt = sort(sample(nrow(data), nrow(data)*.7))
datos_tr<-data[dt,]
datos_te<-data[-dt,]

Función Auxiliar

Esta función toma como argumentos el modelo, los datos de entrenamiento y prueba. Calcula y muestra métricas, y curvas ROC

metricas <- function(modelo, train, test, nombre){
  pred_svm_test= predict(modelo, newdata = test)
  acc_svm_test <- round(measureACC(as.data.frame(pred_svm_test)$truth, as.data.frame(pred_svm_test)$response),3)
  AUC_svm_test <- round(measureAUC(as.data.frame(pred_svm_test)$prob.maligno,as.data.frame(pred_svm_test)$truth,'normal','maligno'),3)
  # ···························································································
  # Predicción TRAIN (naive)
  pred_svm_train = predict(modelo, newdata = train) # por si quiero ver naive sobre training
  acc_svm_train <- round(measureACC(as.data.frame(pred_svm_train)$truth, as.data.frame(pred_svm_train)$response),3)
  AUC_svm_train <- round(measureAUC(as.data.frame(pred_svm_train)$prob.maligno, as.data.frame(pred_svm_train)$truth, 'normal','maligno'),3)
  # ···························································································
  # Cambio el threshold [esto lo hago para train y test]
  acc=NULL
  acc2=NULL
  threshold = seq(0.1,0.95,0.01)
  for (i in 1:length(threshold)) {
          pred = setThreshold(pred_svm_test, threshold = threshold[i])
          acc[i] = measureACC(as.data.frame(pred)$truth, as.data.frame(pred)$response)}
  for (i in 1:length(threshold)) {
          pred2 = setThreshold(pred_svm_train, threshold = threshold[i])
          acc2[i] = measureACC(as.data.frame(pred2)$truth, as.data.frame(pred2)$response)}
  par(mfcol = c(1,2))
  
  new_df1 <- as.data.frame(cbind(threshold,acc))
  new_df1 <- new_df1%>%mutate(sub_data='test')
  new_df2 <- as.data.frame(cbind(threshold,acc2))
  colnames(new_df2) <- c('threshold','acc')
  new_df2 <- new_df2%>%mutate(sub_data='train')
  
  new_df <- as.data.frame(rbind(new_df1,new_df2))
  
  #new_df1[which.max(new_df1$acc),"threshold"] # test 0.39
  #new_df2[which.max(new_df2$acc),"threshold"] # train 0.53
  # ···························································································
  # Gráfico de cómo varía la métrica de performance accuracy, de acuerdo al umbral elegido
  plot_acc <- ggplot(new_df, aes(x=threshold, y=acc)) + geom_line(aes(color = sub_data,linetype=sub_data)) +
          theme +labs(x='Umbral', y='Métrica de performance (accuracy)', 
                       title= 'Evaluación del modelo de Máquinas de soporte vectorial SVM') +
          scale_color_manual(values = c("red", "darkred"),labels=c('prueba','entrenamiento')) +
          scale_linetype_manual(values=c(1,2), labels=c('prueba','entrenamiento')) + 
          labs(color='Conjunto de\n evaluación',linetype='Conjunto de\n evaluación')
 
  print(plot_acc)
  
   # Para independizarnos de la elección del umbral, grafico curvas ROC para las predicciones del modelo SVM con los datos de TEST y TRAIN
  df_svm = generateThreshVsPerfData(list(svm_te = pred_svm_test, svm_tr = pred_svm_train), 
                                    measures = list(fpr, tpr, mmce))
  
  plot_roc <- plotROCCurves(df_svm) + theme +
          labs(title=paste0('Curva ROC del modelo - ', nombre), 
               x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)',
               color='Conjunto de\n evaluación') +
          scale_color_manual(values = c("red", "darkred"), labels=c('prueba','entrenamiento'))
          # geom_label(label="AUC= 0.894", x=0.35, y=0.75, label.size = 0.3, size=4,
          #            color = "red",fill="white") + 
          # geom_label(label="AUC= 0.925", x=0.07, y=0.97, label.size = 0.3, size=4,
          #            color = "darkred",fill="white")
  
  print(plot_roc)
  # ················ Métricas del modelo de SVM ················
  Métrica <- c('valor','datos')
  Accuracy <- c(acc_svm_test,'prueba')
  Accuracy. <- c(acc_svm_train,'entrenamiento')
  AUC_ROC <- c(AUC_svm_test,'prueba')
  AUC_ROC. <- c(AUC_svm_train,'entrenamiento')
  # Imprimo resultados
  kable(rbind(Métrica, Accuracy, Accuracy., AUC_ROC, AUC_ROC.))
  }

Se utilizará la librería mlr, que ofrece herramientas para trabajar con experimentos de machine learning. En el caso de svm, mlr utiliza las funciones del paquete e1071.

Daremos un primer vistazo a un modelo simple con kernel lineal. Para ello se definen la tarea de clasificación y el modelo a ajustar. Posteriormente, se entrena el modelo con los datos previamente reservados par ael entrenamiento.

MÁQUINAS DE SOPORTE VECTORIAL (SVM) con kernel lineal

# Defino modelo SVM
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis") 
lrn_svmL = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost=1)) 
mod_svmL = mlr::train(lrn_svmL, task)

Calcularemos nuestras primeras métricas sobre este modelo, para conocer el desempeño

metricas(mod_svmL, datos_tr[-2], datos_te[-2], 'SVM con Kernel Lineal')

Métrica valor datos
Accuracy 0.826 prueba
Accuracy. 0.872 entrenamiento
AUC_ROC 0.891 prueba
AUC_ROC. 0.924 entrenamiento

¿Cómo se ve la frontera de decisión para este clasificador? Dada la dimensionalidad de nuestro dataset, no es posible observarla de manera directa. A modo ilustrativo, haremos una reducción de la dimensionalidad a través del cómputo de las Componentes Principales.

Visualización de la frontera de decisión para un kernel lineal

# Preparo los datos para análisis de componentes principales
datos_para_acp = data[c(3:7)] # todas las variables numéricas
datos.pc = prcomp(datos_para_acp,scale = TRUE) #escalo los datos
summary(datos.pc)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5
## Standard deviation     1.6043 1.0609 0.7881 0.63208 0.52941
## Proportion of Variance 0.5147 0.2251 0.1242 0.07991 0.05605
## Cumulative Proportion  0.5147 0.7398 0.8640 0.94395 1.00000

Conservaremos sólo las dos primeras componentes, con relativa tranquilidad de capturar con ellas casi el 74% de la varianza de los datos -más que suficiente para un análisis ilustrativo.

data_toplot = datos.pc$x[,c('PC1','PC2')]
data_toplot = data.frame(data_toplot, data[,'diagnosis'])
colnames(data_toplot) = c('PC1','PC2', 'diagnosis')

Definiremos una pequeña función para crear una grilla de datos. Con ella, generaremos datos espaciados uniformemente y cubriendo la totalidad del rango de nuestras variables en la muestra estudiada. Luego, alimentaremos nuestro clasificador con la grilla, y los puntos servirán de testigo para dibujar una frontera de decisión.

make.grid = function(x, n = 75) {
  grange = apply(x, 2, range)
  x1 = seq(from = grange[1,1], to = grange[2,1], length = n)
  x2 = seq(from = grange[1,2], to = grange[2,2], length = n)
  expand.grid(X1 = x1, X2 = x2)
}
set.seed(1)
task = makeClassifTask(data = data_toplot, target = "diagnosis") 
lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost = 1)) 
mod_plot = mlr::train(lrn_plot, task)

xgrid = make.grid(data_toplot[, c('PC1','PC2')])
ygrid = predict(mod_plot, newdata = xgrid)
plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)

# Con esta línea, se pueden conocer los vectores de soporte
#points(data_toplot[,c('PC1','PC2')][mod_plot$learner.model$index,], pch = 5, cex = 2)

En esta figura, se puede ver el diagnóstico real para cada paciente. Y, además, gracias a nuestra grilla podemos ver la clase que el algoritmo asignaría a cualquier otro hipotético punto que quisiéramos incorporar. Los datos están relativamente condensados pero muy superpuestos. Evidentemente, un análisis bidimensional no es suficiente para separarlos -al menos no el obtenido con esta metodología.

Los kernels que ofrece e1071, se pueden conocer a través de la función getParamSet. Además también, si quisiéramos, podríamos investigar los parámetros disponibles para acompañar a cada kernel particular, a través de este mismo objeto.

getParamSet("classif.svm")$pars$kernel$values
## $linear
## [1] "linear"
## 
## $polynomial
## [1] "polynomial"
## 
## $radial
## [1] "radial"
## 
## $sigmoid
## [1] "sigmoid"

¿Cómo se ven las fronteras de separación con otros kernels más flexibles? Si bien no podemos esperar una separación perfecta en dos dimensiones, quizás haya vectores cuyas clases reales sean mejor capturadas por otros kernels.

En primer lugar, observaremos kernels polinomiales de grados entre 2 y 5. Ya puede verse una flexibilidad tremendamente superior al elemental kernel lineal aunque hay algunos puntos difíciles de capturar.

set.seed(1)
for(i in 2:5){
  task = makeClassifTask(data = data_toplot, target = "diagnosis") 
  lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "polynomial", degree=i, cost = 1)) 
  mod_plot = mlr::train(lrn_plot, task)
  
  xgrid = make.grid(data_toplot[, c('PC1','PC2')])
  ygrid = predict(mod_plot, newdata = xgrid)
  plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
  points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
} 
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

¿Cuál es el efecto del parámetro C? Observaremos un barrido entre unos pocos valores del parámetro, para un kernel polinomial de grado 4. Marca la diferencia entre un ajuste más suave o más estricto.

set.seed(1)
for(i in 2:5){
  task = makeClassifTask(data = data_toplot, target = "diagnosis") 
  lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "polynomial", degree=4, cost = i)) 
  mod_plot = mlr::train(lrn_plot, task)
  
  xgrid = make.grid(data_toplot[, c('PC1','PC2')])
  ygrid = predict(mod_plot, newdata = xgrid)
  plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
  points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
} 
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

A continuación, observaremos las fronteras que describe el kernel Radial. Aprovecharemos, simultáneamente, para tener una noción gráfica del efecto del parámetro gamma. Si recordamos la ecuación que define la función del kernel radial, podemos esperar una influencia menor de la distancia entre dos puntos cuanto más grande es gamma

Puede observarse todavía más flexibilidad que en el caso polinomial, y se ve claramente cómo varia la frontera ante el barrido de gamma. Se muestran valores de gamma entre 1 y 6, para un valor constante de C. El efecto de localidad es mayor cada vez.

for(i in 1:6){
  set.seed(1)
  task = makeClassifTask(data = data_toplot, target = "diagnosis") 
  lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "radial", cost=0.8, gamma = i)) 
  mod_plot = mlr::train(lrn_plot, task)
  
  xgrid = make.grid(data_toplot[, c('PC1','PC2')])
  ygrid = predict(mod_plot, newdata = xgrid)
  plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
  points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)
}
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

## Warning in plot.xy(xy, type, ...): NAs introduced by coercion

  #points(datos_tr[,c(x1,x2)][mod_plot$learner.model$index,], pch = 5, cex = 2)

Finalmente, observaremos la frontera del kernel Sigmoide. La curva que dibuja es altamente flexible, pero no ha capturado alunos puntos específicos ¿Será posible obtener un mejor ajuste?

set.seed(1)
task = makeClassifTask(data = data_toplot, target = "diagnosis") 
lrn_plot = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "sigmoid", cost=2, gamma = 0.5)) 
mod_plot = mlr::train(lrn_plot, task)

xgrid = make.grid(data_toplot[, c('PC1','PC2')])
ygrid = predict(mod_plot, newdata = xgrid)
plot(xgrid, col = c("red","green")[as.numeric(unlist(ygrid))], pch = 20, cex = .2)
## Warning in plot.xy(xy, type, ...): NAs introduced by coercion
points(data_toplot[,c('PC1','PC2')], col = as.integer(data_toplot$diagnosis) + 1 , pch = 19)

  #points(datos_tr[,c(x1,x2)][mod_plot$learner.model$index,], pch = 5, cex = 2)

Búsqueda de Hiperparámetros

La definición de C, gamma, el grado -para polinomiales- e incluso el kernel óptimo a utilizar puede hacerse a través de una búsqueda por grilla, aleatoria o bayesiana. A modo de ejemplo, buscaremos el set óptimo de parámetros con una búsqueda aleatoria y cross validation.

El paquete mlr nos sirve en este caso para ejecutar el flujo de trabajo de forma compacta y ordenada. Definida la tarea y el clasificador, definimos el espacio de búsqueda paramétrica. Luego, definimos la estrategia de búsqueda y cv.

# Tarea y clasificador
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis") 
svm <- makeLearner("classif.svm")

# Parámetros a optimizar
kernels <- c("polynomial", "radial", "sigmoid")
svmParamSpace <- makeParamSet(
  makeDiscreteParam("kernel", values = kernels),
  makeIntegerParam("degree", lower = 1, upper = 3),
  makeNumericParam("cost", lower = 0.1, upper = 10),
  makeNumericParam("gamma", lower = 0.1, 10))

# Esquema de búsqueda
randSearch <- makeTuneControlRandom(maxit = 30)
cvForTuning <- makeResampleDesc("Holdout", split = 2/3)

#Parámetros óptimos
tunedSvmPars <- tuneParams("classif.svm", task = task,
                     resampling = cvForTuning,
                     par.set = svmParamSpace,
                     control = randSearch)
## [Tune] Started tuning learner classif.svm for parameter set:
##            Type len Def                    Constr Req Tunable Trafo
## kernel discrete   -   - polynomial,radial,sigmoid   -    TRUE     -
## degree  integer   -   -                    1 to 3   -    TRUE     -
## cost    numeric   -   -                 0.1 to 10   -    TRUE     -
## gamma   numeric   -   -                 0.1 to 10   -    TRUE     -
## With control class: TuneControlRandom
## Imputation value: 1
## [Tune-x] 1: kernel=sigmoid; degree=3; cost=1.27; gamma=0.487
## [Tune-y] 1: mmce.test.mean=0.2471910; time: 0.0 min
## [Tune-x] 2: kernel=polynomial; degree=2; cost=9.45; gamma=5.17
## [Tune-y] 2: mmce.test.mean=0.2359551; time: 0.0 min
## [Tune-x] 3: kernel=radial; degree=3; cost=4.07; gamma=8.97
## [Tune-y] 3: mmce.test.mean=0.1910112; time: 0.0 min
## [Tune-x] 4: kernel=radial; degree=1; cost=1.41; gamma=2.73
## [Tune-y] 4: mmce.test.mean=0.1235955; time: 0.0 min
## [Tune-x] 5: kernel=radial; degree=1; cost=8.45; gamma=1.75
## [Tune-y] 5: mmce.test.mean=0.1573034; time: 0.0 min
## [Tune-x] 6: kernel=sigmoid; degree=2; cost=2.8; gamma=3.5
## [Tune-y] 6: mmce.test.mean=0.3483146; time: 0.0 min
## [Tune-x] 7: kernel=sigmoid; degree=2; cost=3.43; gamma=9.68
## [Tune-y] 7: mmce.test.mean=0.3707865; time: 0.0 min
## [Tune-x] 8: kernel=radial; degree=1; cost=3.64; gamma=5.18
## [Tune-y] 8: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 9: kernel=polynomial; degree=1; cost=2.38; gamma=1.71
## [Tune-y] 9: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 10: kernel=radial; degree=2; cost=1.96; gamma=4.14
## [Tune-y] 10: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 11: kernel=radial; degree=2; cost=9.4; gamma=0.556
## [Tune-y] 11: mmce.test.mean=0.1011236; time: 0.0 min
## [Tune-x] 12: kernel=polynomial; degree=2; cost=9.07; gamma=8.82
## [Tune-y] 12: mmce.test.mean=0.2247191; time: 0.0 min
## [Tune-x] 13: kernel=sigmoid; degree=3; cost=4.09; gamma=5.61
## [Tune-y] 13: mmce.test.mean=0.3370787; time: 0.0 min
## [Tune-x] 14: kernel=radial; degree=1; cost=9.26; gamma=8.8
## [Tune-y] 14: mmce.test.mean=0.1910112; time: 0.0 min
## [Tune-x] 15: kernel=polynomial; degree=3; cost=3.32; gamma=1.69
## [Tune-y] 15: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 16: kernel=polynomial; degree=2; cost=2.29; gamma=0.12
## [Tune-y] 16: mmce.test.mean=0.3820225; time: 0.0 min
## [Tune-x] 17: kernel=sigmoid; degree=3; cost=5.11; gamma=7.29
## [Tune-y] 17: mmce.test.mean=0.3595506; time: 0.0 min
## [Tune-x] 18: kernel=radial; degree=2; cost=3.93; gamma=3.44
## [Tune-y] 18: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 19: kernel=radial; degree=1; cost=6.99; gamma=8.17
## [Tune-y] 19: mmce.test.mean=0.1910112; time: 0.0 min
## [Tune-x] 20: kernel=sigmoid; degree=1; cost=9.07; gamma=4.19
## [Tune-y] 20: mmce.test.mean=0.3483146; time: 0.0 min
## [Tune-x] 21: kernel=polynomial; degree=3; cost=2.49; gamma=3.07
## [Tune-y] 21: mmce.test.mean=0.1685393; time: 0.0 min
## [Tune-x] 22: kernel=polynomial; degree=3; cost=0.237; gamma=5.56
## [Tune-y] 22: mmce.test.mean=0.1685393; time: 0.0 min
## [Tune-x] 23: kernel=radial; degree=3; cost=2.61; gamma=4.17
## [Tune-y] 23: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 24: kernel=radial; degree=1; cost=6.44; gamma=3.53
## [Tune-y] 24: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 25: kernel=radial; degree=3; cost=3.9; gamma=2.36
## [Tune-y] 25: mmce.test.mean=0.1460674; time: 0.0 min
## [Tune-x] 26: kernel=sigmoid; degree=1; cost=6.89; gamma=7.94
## [Tune-y] 26: mmce.test.mean=0.3707865; time: 0.0 min
## [Tune-x] 27: kernel=sigmoid; degree=3; cost=3.12; gamma=7.77
## [Tune-y] 27: mmce.test.mean=0.3707865; time: 0.0 min
## [Tune-x] 28: kernel=polynomial; degree=2; cost=8.35; gamma=7.61
## [Tune-y] 28: mmce.test.mean=0.2247191; time: 0.0 min
## [Tune-x] 29: kernel=polynomial; degree=1; cost=3.37; gamma=8.65
## [Tune-y] 29: mmce.test.mean=0.1573034; time: 0.0 min
## [Tune-x] 30: kernel=sigmoid; degree=2; cost=4.42; gamma=1.85
## [Tune-y] 30: mmce.test.mean=0.3370787; time: 0.0 min
## [Tune] Result: kernel=radial; degree=2; cost=9.4; gamma=0.556 : mmce.test.mean=0.1011236

La salida nos indica la mejor combinación, para esta estrategia y espacio de búsqueda con 30 iteraciones.

[Tune] Result: kernel=radial; degree=1; cost=3.37; gamma=2.96 : mmce.test.mean=0.0786517

A continuación pondremos a prueba el desempeño de varios clasificadores, incluído el que reúne los parámetros óptimos. Contrastaremos su desempeño con el de una regresión logística. Ahora, sobre el problema original y no su versión reducida.

Calcularemos, para cada modelo: 1. Curvas de Accuracy para diferentes thresholds, en sets de entrenamiento y prueba. 2. Curvas ROC en los mismos sets. AUC será nuestra medida de mayor interés.

MÁQUINAS DE SOPORTE VECTORIAL (SVM) con kernel lineal

# Defino modelo SVM
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis") 
lrn_svmL = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "linear", cost=1)) 
mod_svmL = mlr::train(lrn_svmL, task)
metricas(mod_svmL, datos_tr[-2], datos_te[-2], 'SVM con Kernel Lineal')

Métrica valor datos
Accuracy 0.826 prueba
Accuracy. 0.872 entrenamiento
AUC_ROC 0.891 prueba
AUC_ROC. 0.924 entrenamiento

Kernel Lineal: AUC_ROC 0.891 prueba

MÁQUINAS DE SOPORTE VECTORIAL (SVM) con kernel Polinomial

set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis") 
lrn_svmP = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "polynomial", cost=2 , degree=3)) 
mod_svmP = mlr::train(lrn_svmP, task)
metricas(mod_svmP, datos_tr[-2], datos_te[-2], 'SVM con Kernel Polinomial (grado 3)')

Métrica valor datos
Accuracy 0.783 prueba
Accuracy. 0.861 entrenamiento
AUC_ROC 0.895 prueba
AUC_ROC. 0.933 entrenamiento

Kernel Polinomial: AUC_ROC 0.895 prueba

MÁQUINAS DE SOPORTE VECTORIAL (SVM) con kernel radial

Parámetros óptimos: [Tune] Result: kernel=radial; degree=1; cost=3.37; gamma=2.96

# Defino modelo SVM
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis") 
lrn_svmR = makeLearner("classif.svm", predict.type = "prob", par.vals = list( kernel = "radial", cost=3.37, gamma =2.96)) 
mod_svmR = mlr::train(lrn_svmR, task)
metricas(mod_svmR, datos_tr[-2], datos_te[-2], 'Radial')

Métrica valor datos
Accuracy 0.809 prueba
Accuracy. 0.985 entrenamiento
AUC_ROC 0.902 prueba
AUC_ROC. 1 entrenamiento

REGRESIÓN LOGÍSTICA

# chequeo el balance de las distintas clases de  la variable diagnóstico en el conjunto de todos los datos, los datos de prueba y los datos de entramiento.
Entrenamiento <- table(datos_tr$diagnosis) #126 normal, 140 maligno
Prueba <- table(datos_te$diagnosis) # 57 normal, 58 maligno
Total <- table(data$diagnosis) # 183 normal, 198 maligno
kable(rbind(Entrenamiento, Prueba,Total))
normal maligno
Entrenamiento 126 140
Prueba 57 58
Total 183 198 
# Armo modelo de regresión logistica. 
set.seed(1)
task = makeClassifTask(data = datos_tr[-2], target = "diagnosis") 
lrn = makeLearner("classif.logreg", predict.type = "prob")
mod_lr = mlr::train(lrn, task)
metricas(mod_lr, datos_tr[-2], datos_te[-2], 'Logística')

Métrica valor datos
Accuracy 0.817 prueba
Accuracy. 0.872 entrenamiento
AUC_ROC 0.889 prueba
AUC_ROC. 0.925 entrenamiento 
kable(mod_lr$learner.model$coefficients)
x
(Intercept) -5.4691572
edad 0.0692457
creatinina -1.0962220
LYVE1 0.4330571
REG1B 0.0027155
TFF1 0.0020305

COMPARACIÓN DE MÉTODOS DE CLASIFICACIÓN SUPERVISADA

Finalmente, compararemos las métricas calculadas sobre el desempeño de todos los modelos construídos.

# ······················ curvas ROC para todos los modelos ········································
pred_lr = predict(mod_lr, newdata = datos_te[-2])
pred_svmL = predict(mod_svmL, newdata = datos_te[-2])
pred_svmP = predict(mod_svmP, newdata = datos_te[-2])
pred_svmR = predict(mod_svmR, newdata = datos_te[-2])

df_todos = generateThreshVsPerfData(list(lg=pred_lr, svmL = pred_svmL, svmP = pred_svmP, svmR = pred_svmR), measures = list(fpr, tpr, mmce))

plotROCCurves(df_todos) + theme + 
        labs(title='Curvas ROC de modelos de clasificación supervisada (datos de prueba)',
             x='Tasa de falsos positivos (FPR)', y='Tasa de positivos verdaderos (TPR)', 
             color=' Modelo en\n evaluación') +
        scale_color_manual(values = c("red", "black", "blue", "darkgreen"),
                           labels=c('Reg log','SVM (lineal)','SVM (polinomial)','SVM (radial)'))+
        theme(legend.position=c(0.915,0.25))

# ············ valores AUC para todos los modelos cuando se consideran todas las variables ···················
AUC_lg_te <- round(measureAUC(as.data.frame(pred_lr)$prob.maligno, as.data.frame(pred_lr)$truth, 'normal','maligno'),3)
AUC_svm_teL <- round(measureAUC(as.data.frame(pred_svmL)$prob.maligno, as.data.frame(pred_svmL)$truth, 'normal','maligno'),3)
AUC_svm_teP <- round(measureAUC(as.data.frame(pred_svmP)$prob.maligno, as.data.frame(pred_svmP)$truth, 'normal','maligno'),3)
AUC_svm_teR <- round(measureAUC(as.data.frame(pred_svmR)$prob.maligno, as.data.frame(pred_svmR)$truth, 'normal','maligno'),3)

AUC_values <- rbind(AUC_lg_te, AUC_svm_teL,  AUC_svm_teP,  AUC_svm_teR)
AUC_values <- as.data.frame(AUC_values)
AUC_values$Modelo <- c('Reg Log','SVM L','SVM P','SVM R')
colnames(AUC_values) <- c('Area debajo de la curva (AUC)','Modelo')
row.names(AUC_values) <- NULL
AUC_values <- AUC_values%>%dplyr::select(2,1)
# Imprimo resultados
kable(AUC_values)
Modelo Area debajo de la curva (AUC)
Reg Log 0.889
SVM L 0.891
SVM P 0.895
SVM R 0.902

Con un AUC de 0.902 sobre el set de test, el Kernel radial es superior. Las diferencias en rendimientos, con esta configuración, no son demasiado grandes. El paso siguiente, sería comparar el rendimiento luego de una búsqueda exhaustiva de hiperparámetros para todos los modelos.

Específicamente, el “costo” principal de SVM tiene que ver con la pérdida de interpretabilidad directa.

Nota Bibliográfica Snippets de código fueron basados en: 1. https://www.r-bloggers.com/2019/10/support-vector-machines-with-the-mlr-package/ 2. https://www.datacamp.com/community/tutorials/support-vector-machines-r